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We explore the localization of a quasi-one-, quasi-two-, and three-dimensional ultra-cold gas by 
a finite-range defect along the corresponding 'free'-direction/s. The time-independent non-linear 
Schrodinger equation that describes a Bose-Einstein condensate was used to calculate the maximum 
non-linear coupling constant, Qmax, and thus the maximum number of atoms, Nmax, that the 
defect potential can localize. An analytical model, based on the Thomas-Fermi approximation, is 
introduced for the wavefunction. We show that Qmax becomes a function of RoVVo for various one-, 
two-, and three-dimensional defect shapes with depths Vb and characteristic lengths 7?o- Our explicit 
calculations show surprising agreement with this crude model over a wide range of Vo and _Ro- A 
scaling rule is also found for the wavefunction for the ground state at the threshold at which the 
localized states approach delocalization. The implication is that two defects with the same product 
RoV^ will thus be related to each other with the same Qmax and will have the same (reduced) 
density profile in the free-direction/s. 

PACS numbers: 03.75.-b, 03.75.Hh, 05.30.Jp, 67.85.-d, 67.85.Bc 



I. INTRODUCTION 

Since the experimental realization of Bose-Einstein 
condensates (BEC) in 1995 interest in the physics 

of ultracold atoms has grown and new areas of research 
have emerged. At the same time, it has renewed the in- 
terest in studying the collective dynamics of macroscopic 
ensembles of atoms occupying the same single-particle 
quantum state [§-[1]. This, in turn, has created the need 
for new technology to study ultracold atoms. There are 
two technologies that motivate our present study; one is 
the atom chip, the second is sculptured optical field-based 
atom trapping. Both of which are capable of generating 
a myriad of trap geometries. 

Our fundamental goal in this work is to determine the 
general scaling law that determines the maximum num- 
ber of atoms that can be trapped by attractive defects de- 
scribed by 1-D, 2-D, and 3-D potentials when combined 
with confining harmonic atom traps. This corresponds 
to determining the limit at which the system transitions 
from a bound to a scattering state. 

Our first motivation are the atom chips which play an 
important role in atomic physics, enabling the cooling 
and trapping of a BEC in a waveguide which is created 
by magnetic fields generated above patterned micro- wire 
circuits @. Atom chips have already enabled the study 



of matter- wave interference phenomena Q , and may im- 
prove other atomic measurement devices such as atomic 
clocks in the future. Ideally, the BEC is loaded/trapped 
in the transverse ground state, but is allowed free-space 
propagation along the third dimension. An ultra-cold 
wave-packet can then be transported through quasi- 1- 
D waveguides due to the strong confinement in the two 
transverse dimensions Q. 

Our second motivation is the experimental realization 
of potential traps with different shapes that have been 
achieved using a rapidly moving laser beam that paints 
a time-averaged optical dipole potential. There, a BEC 
is created in arbitrary geometries 0, ■ Effectively, the 
BEC confinement can be strong in one-dimension, re- 
ducing the BEC to be quasi-2-D in shape. The BEC 
can then be manipulated in these two-dimensions using 
a time-averaged potential. Both the atom chip and sculp- 
tured traps motivate us to explore a variety of shapes for 
a basic trap which includes a defect potential. 

In the study of the structure of a gas of ultra-cold 
atoms, it is necessary to account for the interaction of 
the N atoms. That is, each atom moves in an average 
field due to the other — 1 atoms that surrounds it. 
The properties of a BEC at T = is well-described by 
a mean-field approximation which results in a non-linear 
Schrodinger equation (NLSE) for the single-particle or- 
bitals lUi] 
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The nonhnear couphng constant characterizes the short- 
range pairwise interactions, and is given for N bosons 
by 93 = 47ras(iV — l)/a± according to (number con- 
serving) Hartree-Fock theory [ll|. This constant de- 
pends on the s-wave atom-atom scattering length, a^, 
of two interacting bosons. Note that alternative treat- 
ments give 33 oc N, and then the above NLSE is known 
as the Gross-Pitaevskii equation (GPE) Q. Eq. ^ 
has been written in oscillator units (o.u.), with lengths 
a_L = ^hjmujx^ for a chosen frequency wj^ and m being 
the mass of the individual atoms that compose the BEG. 
The energy is thus given in units of the oscillator en- 
ergy ?ia;_L, and time in units of Back in S.I. units, 
gf = ^-Kh^asiN - l)/m = h^gi/iai^m). 

For the present study, we chose time-independent po- 
tentials for V{y) with one of three structures: 

Vi{x, y, z) = imc^i(y2 + -J") + Vi{x) , (2) 

V2{x,y,z) = ]^muj^{z^) + V2{x,y) , (3) 

V3{x,y,z) ^Vl{x,y,z) . (4) 

That is, we partition the defect potential V' , which is of 
interest here, from the transverse harmonic trap potential 
applied at a frequency wj,. 

The VI defect potentials, for example, can be gener- 
ated by a local modification of the transverse waveguide 
confinement, such as a constriction or a local cur- 

vature |16l - [l8| . Whether a defect potential acts as an 
obstacle or a sink, a wave-packet will interfere with, and 
possibly lose atoms as it goes through the defect due to 
the non-linearity [iqI . [2oI |. changing the interaction for 
any subsequent atom. In general, propagation through 
a perturbation in a 1-D waveguide results in unwanted 
transverse excitations of the BEG 0, [2l| . 

The presence of a transverse-wj, potential allows for 
the reduction of the NLSE from 3-D to either a 2-D or 
1-D form. In the 1-D, or Vi, case we assume that the 
tight-waveguide limit applies where the longitudinal- (a;) 
size of the BEG is larger than its transverse- (y, z) cross 
section. This allows for us to integrate out the trans- 
verse dimensions which are energetically frozen in the 
harmonic ground state. This results in a 1-D NLSE, 



19^ 



+ V[{x)+gim' 



^{x,t) (5) 
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with an effective 1-D coupling constant, gi 
2as{N — l)/a±. Back in S.I. units gf^ — g^' 
and note that the intermediate reduction from 3-D to 
2-D gives a similar NLSE with constant g2 = 33/ V2n. 

By systematically adding more and more atoms with 
tts > into the potential, at some gs^max the total 3-D 
eigenenergy, £3, of the system becomes null with respect 
to the transverse trap energy [l^. For the three geome- 
tries in Eq. this corresponds to £3 = fiuj±, £3 — ^fku^ 
or £3 = 0. Determining the ga^max is non-trivial since at 
this point the system is no longer bound, the delocalized 



atoms can reach |r| = 00, and thus cannot be repre- 
sented by a square-integrable wavefunction, and needs to 
be considered as the limit £3 — >■ from below. 

The first analytical approximations for gi^max were si- 
multaneously derived by Garr et al. |22| and Leboeuf and 
Pavloff [3. Garr et al. solved the 1-D NLSE for a fi- 
nite square-well finding localized solutions when the 1-D 
eigenenergy, Si < 0. The transition at £1 = was found 
in terms of an approximate expression for gi^max 

contain- 
ing i?o\A^ terms where the potential well has depth, Vq, 
and width, 2Ro. (see Eq. (21) of Ref. jUl) which our 
results will show agreement with. 

Leboeuf and Pavloff [l6l | derived approximate expres- 
sions in terms of the maximum number of atoms that a 
1-D potential could support based on the area enclosed 
by the potential, A = | V{{x)dx\. In the low-density 
BEG limit, this translates from their h = m = 1 units to 
gi,max = 2A when in o.u.. They argued that, in general, 
this was expected 'to be very accurate' despite using a 
V({x) — >• —XS{x) approximate mapping. They did not 
explicitly validate their formula against explicit calcu- 
lations. More recent analytical work by Seaman et al. 
[HI showed that, for a potential V({x) — ~l3S{x), the 
gi.max = 4/3 exactly. Our results agree with Seaman et 
al., in that, there is a factor of 2 underestimate in the 
treatment of Leboeuf and Pavloff in the limit of weak 
potentials. Our results also agree with Garr et al. [l^l 
which show that the Leboeuf and Pavloff expression has 
a factor of 2 overestimate in the limit of strong potentials. 

The final paper of direct relevance to the present study 
is that of Adhikari Q who solved the 3-D NLSE with 
a finite spherical-well potential, and computed the max- 
imum nonlinear coupling constant, gs^max- 

There, a 3-D 

variational ansatz was applied and an almost linear rela- 
tion was plotted (see their Fig. (3)) in which gs^max oc Vq, 
the depth of the well. Adhikari noted that the variational 
ansatz underestimates the magnitude of g^^max, but gave 
no calculations to quantify this. 

Furthermore, of particular interest is that Adhikari [23| 
also applied a Thomas-Fermi approximation (TFA) to 
compute, in our notation, g^^max = 47ri?oVb/3. In their 
paper it is noted the TFA "is inadequate for calculat- 
ing" 53, max, leading to a "much smaller" value than the 
variational ansatz. It was thus surprising to us when we 
applied the same TFA in preliminary 1-D calculations 
25], and found excellent agreement with explicit calcu- 
lations of gi^max- In the present paper we find that the 
TFA gives g2,max and gs^max's accurately for all of the 
potential defects that we consider. Our 3-D results show 
that both the gs^max results and conclusions of Adhikari's 
variational ansatz [2^ are not accurate. 

In this work, we present some analytical results in Sec. 
|TT] using the TFA for i = 1-, 2-, and 3-D defect potentials 
and deduce universal scaling rule properties for gi,max 
and the wavefunction (density profiles). The gi^max is 
seen to follow a scaling law that depends on Rq^/Vq, 
across a wide range of width parameters and potential 
shapes. In the 1-D case, we present results for a square- 



well, triangle, truncated harmonic, truncated double har- 
monic, and truncated half circle defect. We have previ- 
ously published results of four of these trap shapes solely 
in 1-D for gi^max jssj - Here, we expand and extend that 
work into 2-D and then 3-D for a wider range of potential 
defects. For the 2-D traps case, we study a rectangle well, 
truncated harmonic and pyramid defect. For a 3-D trap 
we use a spherical, rectangular cuboid and a truncated 
harmonic well defect. Furthermore, we report a scaling 
rule for the wavefunction in terms of the range of the po- 
tential. In Sec. mil we present our numerical approach to 
solve the NLSE based on a finite-difference method on a 
numerical lattice and a Gauss-Seidel procedure to obtain 
the numerical ground state on the lattice for the 1-, 2-, 
and 3-D NLSE. In Sec. llVi we present and discuss our 
results. Finally, in Sec. |V]we present our conclusions. 



II. ANALYTIC RESULTS USING THE 
THOMAS-FERMI APPROXIMATION 

The TFA considers the limit of strong interactions be- 
tween atoms that form an ultra-cold BEG and allows for 
some useful expressions for the single-particle wavefunc- 
tion to be obtained [1^, [27| . A BEG is said to be in the 
Thomas-Fermi (TF) regime when the interaction energy 
dominates over the zero-point energy [23] . The TF 
states are strictly localized by a potential, and thus their 
behavior was not initially expected to be able to mimic 
the NLSE localized states as they approach delocaliza- 
tion as was previously noted by Adhikari Instead, 
we will show that this approximation gives us a general 
(JJoaA^) -based scaling law that agrees with the numeri- 
cal NLSE solutions. 

The TFA neglects the kinetic energy in the NLSE, and 
therefore the time- independent NLSE in either is 1,2,3 
dimensions becomes 

[Vl{v) + gJ^\^TF{v)\^] *TF(r) = e,-^TF{v), (6) 

where Si is the z-D single-particle orbital energy. The 
maximum nonlinear coupling constant due to the defect 
potential occurs when — 0, so Eq. ^ reduces to 



^'TF(r) = lim 



TF 
9i,max 



(7) 



Through normalization, |5'(r)p(ir* — 1, one obtains 
that gf,max satisfies the relation 



Vl{v)dr' 



(8) 



where Q. is the contour around the defect where y/(r) = 
0. Thus, the non-linear coupling parameter gf^^^^ just 
depends on the extension fl (area or volume for the 2- or 
3-D case) of the defect potential in the TFA. Hence, for a 
given scattering length a^, atom mass m and transverse 
frequency lo±, the corresponding maximum number of 
atoms trapped by a defect, Nmax, can be determined. 



A. One-dimensional case 

We consider here the following five 1-D defect potential 
shapes. In all of these cases, Vb > is the strength of the 
potential and 2Rq is the width of the potential region. A 
square potential 



^^'^^"^'-X ,\x\>Ro, 



(9) 



a half circle potential trap 



1^1 < ^0' (10) 
,|2:|>i?o. 

a (truncated) harmonic potential trap 



-Vn 



1 



\x\ < i?n 



10, 

\x\ > Ro, 

a symmetric double harmonic potential trap 



4^0 Ul +1^ 



VIM = I 







, -Ro < X <0, 

,0 < X < Rq. 
, \x\ > Ro, 



(12) 



and a triangle potential 

-^0 ( 1 + ^ 



-^0 ( 1 - 




, -i?o < X <0, 
,0 <x < Ro, 
, \x\ > Ro, 



(13) 



From Eq. (|5]), each potential gives, respectively, a non- 
linear coupling constant given by 



TF 



2 ( -Ro Vn) ) /Rq, square — well 

n{RoVVo?/i2Ro): circle 

4(i?o\/K))V(3^o), harmonic 

4(i?o\/K))V(3-Ro), double 

[RoVVof/Ro, triangle. 



(14) 



Thus, it seems that gf^ax has a universal dependence on 
RqVV) in the 1-D case, independent of the shape of the 
defect potential, i.e gi = f{RoV^)- 

A more general, but still approximate, expression for 
the nonlinear coupling constant that includes the kinetic 
energy term in the NLSE has been found by Garr et 
al. for the one-dimensional square- well, (see Eq. (21) 
in Ref. \2^), and is given by 



51, 



.(7)-f 

Ho 



%/2 



-\/27 



1 + 2e-v^'/ - e-2\/27 



-27 



(15) 

where 7 — RoVVq. In the limit when 7^1, the 
exponentials can be neglected and Eq. reduces 
to gi^maxin) = 2(i?oV^)^/-Ro. This agrees with the 



TFA found for the case of the square-well (first line in 
Eq. (HH)). The other defect potentials considered here 
appear to have no equivalent analytical expression in the 
literature. 



B. Scaling of the 1-D NLSE 

Due to the chosen geometries of the defects, one can 
make an additional change of variables to x = x/Rq and 
El = ei/Vo obtaining the following reduced 1-D NLSE, 



2VoRl 9x^ 



v{x) + ^m' m = eMx), 

voRq 

(16) 

where $(a;) = •\/i?o5'(x) such that $(a;) remains nor- 
malized in the x range. Here, V{x) = V({x)/Vo is the 
reduced defect potential with a maximum strength of —1 
and defined only in the range \x\ < 1. 

For large values of (i?oV^)^i one notices that the ki- 



netic energy term in Eq. ([T6| can be neglected, justifying 
the TFA when Rq and/or Vq are large. The non-linear 
term can not be neglected since gi = f{Ro\/\^), as pre- 
viously mentioned. The equivalent of Eq. ([8]) in these 
reduced units can be thought of as a shape factor, a, 
given by the area under the reduced potential: 



aiD 



9i 



TF 



(i?oVT^)2 



V{x)dx 



(17) 



|S|<1 



This shape factor takes the values of 2, 7r/2, 4/3, 4/3, 
and 1 for each 1-D potential, respectively. 

On the other hand, from Eq. p^ . we now note that 
defects of the same type with the same Rq\/Vq have the 
same gi^max and the same reduced wavefunction, there- 
fore, there is a fundamental relation between the density 
profile for different parameters. We will discuss below 
some further examples where this scaling law holds. 



a truncated harmonic (parabolic) well with circular base 
^^^^^y)'\ , otherwise, 

(21) 

thirdly, a pyramid (triangle) potential 



Vi{x,v) = I 



-Vo 



iV+Ry) 



,y<\x\, -Ry <y<0 
,y>\x\, Ry>y>0 
,\y\<x, -R^ <x<Q (22) 



Rn 

(y~K) 
Ro 



T . {x~R^) 
^0—R^ 





\y\> X, R^> x>0 
, otherwise 



where each line represents the equation on each pyramid 
wall. 



D. Scaling of the 2-D NLSE 

Following a procedure similar to the 1-D case, the 2- 
D NLSE can be further rewritten in terms of the re- 
duced variables. Making the change of variables: x = 
x/y^R^Ry, y = y/y^R^Ry, $ = ^R^Ry"^, and £2 = 
S2/V0, one obtains 



1 



2VoRxRy 



dx^ dy"^ 

92 
VoRxRy 



+ V^{x,y)+ (23) 
|$p}$(fc,y) = e2<^{x,y) 



where V2 = is the 2-D reduced potential. For 

example, for the rectangular well, we have 



V2{x,y) 



-1 ,if \x\ < 




and \y\ < 



, otherwise . 



(24) 



C. Two-dimensional case 

For the case of a two-dimensional traps, the time 
independent NLSE is given by 

dx^ dy'^ 

g2\-i'\^}'i>{x,y)^e2'if{x,y). 

Thus, the application of the TFA again gives the follow- 
ing expression for the 2-D coupling constant at delocal- 
ization 



+ V^{x,y) + 



(18) 



TF 

c)2,max 



V2{x,y)dxdy. 



(19) 



In this case we have considered the following three 2-D 
defect potentials: a rectangular well 



V^{x,y) = 



-Vq , if \x\ < R^ a.nd \y\ < Ry 
, if |x| > Rx and \y\ > Ry , 



(20) 



This gives the following scaling law and expressions for 
the nonlinear 2-D coupling constant in terms of the 2-D 
shape factor 



a2D 



TF 
92 

VaRxRy 



V2{x,y)dxdy (25) 



where the limits on integration are for \x\ < Rx/ Ry 
and \y\ < yjRy/Rx- For the defects considered here. 



TF 

J2,max 



AVoRxRy , Rectangle 
§VqRI , Parabohc 
I Vq Rx Ry , Pyramid . 



(26) 



such that the three shape factors are 4, tt/2, and 4/3 
and thus g2 — fiVoRxRy) which for a symmetric case 
becomes 52 = f{RoV^)- Again, large values of V^RxRy 
can justify the TFA since the kinetic energy term can be 
neglected in Eq. (j23p . but not the non-linear term. 



E. Three-dimensional case 

For the case of three-dimensional traps, the time- 
independent NLSE in oscillator units is given by 



{-I 



,2 n 



+ Vi{x,y,z)+ (27) 



&^ 

53|^'|2|^'(a;,y,z) =e3*(x,2/,z). 

In this case we have considered three different 3- 
dimensional defect potentials. Firstly, a spherical well. 



-Fo ,if r < i?o 
, if r > i?o 



where r — -^/x^ + j/^ + z^. 
cuboid well. 



(28) 

Secondly, a rectangular 



yi{x-,y,z) 



-K) ,if \x\ < R., \y\ < Ry, 

and \z\ < Rz 
,if |a;| > i?^, \y\ < Ry, 
and l^l > Rz , 



(29) 



(which forms a cube when R^ = Ry — Rz). Thirdly, a 
truncated harmonic (parabolic) well, 

-K)(l-rVi?2) ,ifr<i?o 







, otherwise , 



(30) 



where r = \J -\- y"^ -\- z^ 

The application of the TFA again gives the expression 
for the three-dimensional constant, 



V^'(x, 7/, z) dxdydz . 



(31) 



F. Scaling of the 3-D NLSE 

Repeating the procedure used in the 1- and 2-D case, 
the 3-D NLSE can be rewritten in terms of reduced 
variables. Making the changes, x — x/{RxRyRzY^^ ^ 
y = y/iR^RyRzY^^ and z = z/{RJiyRzY'^ for the 
reduced coordinates and $ — RxRyRz"^ for the re- 
duced wavefunction and £3 = es/Vo for the energy, one 
obtains 

dx"^ dy"^ dz^ 



1 



2Vo{RxRyR.)^/^ 
+V3{x,y,z) + 



(32) 



gi/iRxRyRzY'^ 



x,y,z) 



VoiRxRyRzf/^ 

= £3^{x,y,z) 

where V3 = Vj'/Vb is the 3-D reduced defect potential. 
For example, for the rectangular cuboid well, we have 

/ jj2 \ 1/3 

%{x,y,z)={ l^l<(iSr) ' (33) 

aild 1^1 < (-Rfft;) 

, otherwise . 



1/3 



In this case, the TFA gives the following scaling law 
and expressions for the nonlinear 3-D coupling constant 
in terms of the 3-D shape factor 



TF 

93 ,max 

VoRxRyRz 



V3{x,y,z)dxdydy , (34) 



where the integration is over the range \x\ < 
{Rl/RyRzY/^ \y\ < {Rl/RxRzY'\ and \z\ < 

{Rl/RxRy)'/\ 

For the case of a spherical well, a rectangular cuboid, 
and harmonic parabolic well, we have that the TFA gives 



TF 



^^V^Rl/i 

93,'max = { SVoRxRyRz 

'^VoRl 



, Spherical 
, Cuboid 
, Parabolic , 



(35) 



such that the shape factor is 47r/3, 8 and 87r/15, 
respectively. Similarly to the 1- and 2-D case, 
.93 = f (VoRxRyRz), which for a symmetric case 
g3,rnax/Ro = .f{RoV^), such that whcn large values of 
Vo{RxRyRz)^^^ are involved, then the TFA is justified 
in Eq. ([5^ . Remember, the non- linear term can not be 
neglected. 

Thus, in summary, defining the defect shape factor 



(36) 



where i is the dimension of the space and the integra- 
tion is on the ft space that contains the trap. Then the 
maximum nonlinear coupling constant that a defect can 
support is given by 



ctiVoRx 



i = ID 
i = 2D 
UiVoRxRyRz, « = 3D , 



li.^.ax = { OiiVoRxRy 



(37) 



and, for the particular case of symmetric traps where 

Rx = Ru = Rz = -Rn, then 



^i.max 



a^iVnRoy^Ro, J = 3D. 



That is. 



a,{V%RQ)yRo, i = lB 
a.iVVoRo?, « = 2D (38) 



(39) 



fiRoVVo) 



for the i-dimension case and the larger the shape factor, 
the larger the number of trapped atoms by the well. 

In order to verify these expressions outside of the TFA, 
let us solve Eqs. (HH), and (|32]) for the reduced 

wavefunction, by a numerical procedure. In the next 
section we will outline the numerical method we im- 
plemented to solve the NLSE for one, two, and three- 
dimensions for solutions approaching the point of wave- 
function delocalization, that is, when — ^ 0. 



III. NUMERICAL APPROACH 

In this section we present the two computational meth- 
ods we used to compute the wavefunction, and the algo- 
rithm we used to determine the Qmax- 

A. Crank-Nicolson method 

By using finite-differences and the Crank-Nicolson 
method (CN) [28l - l30| . the ground state solutions to the 
time-independent NLSE can be found for a given non- 
linear coupling constant g without approximation. To 
do so, the time-dependent NLSE is evolved in negative 
imaginary time [3l|. The kinetic and potential energy 
terms can be efficiently evolved by means of a symmetric 
split-operator method [11] 

Here T is the kinetic energy operator and V is the po- 
tential energy operator in the NLSE. 

This requires that the wavefunction to be discretized 
in space on a numerical grid, viz. $(5;^, j/j, Zk,tn) ^i)k- 
The 1-D NLSE in this approach, for example, becomes 

{^k - i^{^k+i - 2£_k + £.k-i)} = 

{/fe + K,/fe+i-2/fc + /fe-i)}, (41) 

where v — «At/(4Ax^Voi?o)i ^^'^ PO" 
tentials are in fk — exp{— ?AtV^"/2}$jJ and 

a = exp{iAtF,"+V2}$r' with T4" = Vk + 
5ii?o/(Vbi?g)|$fcp. Eq. gj) can be written in 

matrix form as A+^"+^ = A^^. Note that is a 
constant matrix for fixed At and Aa;. For the multi- 
dimensional problem, the unitary operators are applied 
in sequence, e.g. e^'^*"^ — e^*'^*'^='e~*'^*'^ye~*^*'^^, 
thus only requiring (many) tri-diagonal matrix solves at 
intermediate stages. 

B. Gauss-Seidel method 

We also implemented the Gauss-Seidel (GS) method 
[3^ in 1-D, 2-D and 3-D to find the ground state solution 
of the time- independent NLSE. In this case, the energy 
is evaluated given an improved solution at each point of 
the numerical grid. The GS is much simpler than the 
CN method, and serves as a valuable cross-check that 
our solutions are converged. 

For example, in the 1-D case, the wavefunction evalu- 
ated in the k-th grid point, $fc, is replaced by where 

*;«(!- /?)$,. + (42) 

/3(^fc+i + $fc-i) 

2[1 + {VkVo + g^^/Ro - e,)(i?oAx)2)] ' 

where f3 is the relaxation parameter that ensures conver- 
gence to the lowest energy state. In our case we take 




X (o.u) 

FIG. 1. (Color on-line). The computed NLSE ground state 
wavefunction for a 1-D square-well defect. The {+) symbols 
are the results of the CN method and (A) for the GS method 
in a square-well with gi = 25, Vo ~ 50, and Ro = 1.0 o.u. 
which gives ei = —35.41668 o.u.. The analytical wavefunction 
solution (solid line) of Ref '22\ , and the TFA (dashed line) 
of Eq. ((T]) are also shown for comparison. 

/3 = 3/4 as a compromise between convergence time and 
precision [s^. Extension of the GS method to 2-D and 
3-D is straightforward. 

The results of both the GS and CN computational 
methods are illustrated in Fig. [TJ This shows the wave- 
function, ^'(x), for a 1-D NLSE with a strong non- 
linearity {gi = 25 and £1 = —35.41668), that is trapped 
by a strong square-well (Vb = 50 and Rq = 1.0 o.u.). 
Note that these are all given in oscillator units, and not 
the reduced units. This is compared with the analytical 
wavefunction of Carr et al. (see Eq. (12) in Ref. ,22|). 
The excellent agreement validates our numerical results. 
Also shown is the TFA of Eq. ([7]), where the wavefunc- 
tion takes the same shape as the defect and is identically 
zero outside the well. 



C. Determination of gi^max 

The main computational challenge is to systematically 
increase g to find gi^max for a given Vq and Rq for i = 1- 
, 2-, and 3-D potentials. Thus, the problem is reduced 
to a root search for which ei{gi^max) = for a given Rq 
and Vq. As gi is increased, the wavefunction penetrates 
further into the classically forbidden region, until gi^max 
is reached, at that point the solution determined is no 
longer a localized wavefunction [2i| . To determine gi,max 
(and hence Nmax) numerically, it is required that the dis- 
cretized grid be large enough to contain the very slow de- 
cay of the wavefunction when gi approaches gi^max- Note, 
however, that we strictly use ^{±Xmax) = 0, where we ef- 



fectively have a box Vi{±Xmax) — ^ oo. This means that 
our Ei > scattering states are artificially constrained 
and are always square integrable. 
To ensure that the determined cji^rnax 

is accurate to the 

precision that we demand, the results simply need to be 
insensitive to the location of the boundary. In our 1-D 
cases we could simply choose Xmax = 300 and Aa; = 0.01. 
The 9i,max^s were determined by choosing the defect pa- 
rameters Vq ranging from 0.01, 0.1, 0.5, 1.0, 5.0, 10.0, 
20, and 50.0 o.u., while for the width we choose Rq from 
0.5, 1.0, 2.0, 5.0, 10.0, and 20.0 o.u. commensurate with 
physical parameters of typical BEC experiments (as dis- 
cussed later in Sec. IIV E[) . 

For the 2-D cases we chose Xmax — Vmax = 5.0 
in the reduced units space where the potential reaches 
only unity range and Ax — Ay — 0.05 which gives a 
200x200 grid points in the wavefunction. For the 3-D 
cases, we similarly used x^ax = Vmax = Zmax = 5.0 
and Ax = Ay = Az = 0.05 with a wavefunction size of 
200 X 200 X 200. Thus, increasing the dimension in the 
calculation introduces a large memory footprint, as well 
as many tri-diagonal matrices to solve in the CN calcula- 
tions. Compared to the 1-D calculations, we restrict the 
range of defect parameters explored to ensure accuracy, 
whilst still spanning the parameter regimes of interest. 

In all cases, we assumed convergence for a given gi 
when Si from one iteration to the next changes less than 
Asi = 10^^. The gi was then incremented until Si — 
was located to within lO"'^. The amount of increment of 
gi was 1% of the TFA initial value. When > was 
obtained for a given gi, a step back was performed and 
Agi was reduced in half, until we reached = from 
below. 



IV. RESULTS AND SCALING LAWS 
A. The 1-D square well case 

The numerical and analytical results for the 1-D square 
well are shown in Table H] for the smallest and largest 
value of Rq used in this work. The gi^max value was 
firstly determined for a range of defects spanning large 
and small values of the product Roy/\^. We have chosen 
to translate this into the maximum number of trapped 
atoms that this corresponds to, specifically for a cloud 
of ^^Rb atoms (see caption). These are then compared 
against the expressions from the TFA, and those obtained 
by Carr et al. [22*1 and Leboeuf and Pavloff Tg^ . 

For large Rq^/Vq, the potential is strongly binding and 
our calculations agree closely with both the TFA and 
that of Carr et al. |22|. In this limit, as we increase 
the gi towards the gi^max value, the numerical wavefunc- 
tion resembles the TF wavefunction until we approach ex- 
tremely close to the gi^max and the wavefunction rapidly 
delocalizes. Thus, our numerical method for computing 
gi,max is quite accurate in this regime and not affected 
by the presence of the boundary conditions. The formula 



TABLE I. The maximum nonlinear coupling constant, g-i^max, 
for 1-D square well potentials with various widths (27?o) and 
depths (Vo) as determined by the 1-D numerical calculations. 
The corresponding maximum number of atoms trapped by 
the potentials are given as the Ni^max columns assuming that 
the atomic species is *'^Rb (with as = lOOao, where ao is 
the Bohr radius) trapped by a transverse frequency of cj = 
2tt X 100 rad/s (thus a± — 2.7 microns). The superscripts 
denote (a) the numerical calculation, and (b) the present TFA 
from Eq. (1141) . The approximations of (c) Carr et al. j22i] and 
(d) Leboeuf and Pavloff [3 are given for comparison. 
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0.5 


0.05 


0.044346 


12.321 


13.764 


48.597 


26.528 


0.5 


0.1 


0.141031 
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52.056 


0.5 
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0.852832 
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0.5 
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10 


0.5 
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2553.8 
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10 


1 


21.023785 
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5106.6 


10212 


10 
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102.877445 


26264 


25529 


25529 


51057 


10 


10 


204.679689 


52252 


51057 


51057 


102113 



derived by Leboeuf and Pavloff [16[, with weak binding 
potentials in mind, has a factor of 2 overestimate in the 
strong binding limit. 

It is also interesting, however, to reconcile the values 
in Table U in the small Rq\/Vq limit. Seaman et al. (23j 
gives a factor of gi^max = 4A, the same as the Carr et 
al. [23 in the limit of a small perturbation. That is, 
the system should approach that of binding to a delta- 
function of 'strength' A — 2RqVo. We expect that our 
-Ro ~ 0.5, Vo = 0.05 numerical calculation should be close 
to gi^max = 4 X 2RoVq = 0.2 and not the gi^max = 0.05 
which the TFA estimates. For this weakly bound case we 
find that, even with our numerical wall located at Xmax = 
±300, our gi^max and thus Ni,max are underestimated 
for the weakest traps. Essentially, the overall energy of 
the system is artificially raised due to the system being 
trapped in R finite sized grid, and thus the Qi^rnax 

needed 

to reach delocalization is significantly underestimated for 
small Rq\/Vq. Finally, Leboeuf and Pavloff [l^ give a 
factor of 2 underestimate for small perturbations even 
though they assumed the limit of a small perturbation in 
their derivation. 



B. The general 1-D trap cases 

The results for gi^max for all five 1-D defects are shown 
in Fig. [5] as a function of Rq^/\^. In order to avoid over- 
lapping, the results for the different shape defects have 
been scaled by a factor of 10 between each other in the 
ascending order as given by Eq. (IT4l) . In the same figure 
we show the Thomas-Fermi approximations represented 
by the solid lines, while the long-dashed line is a guide to 
the eye that follows the numerical results. For the case 




Ro(Vo)"' 

FIG. 2. 1-D scaling law for gi,max weighted by Ro for a square 
(a), circular (b), harmonic (c), double harmonic (d), and tri- 
angular well (e) defect, respectively, of various lengths Ro and 
depths Vq. The numerical results are given as various sym- 
bols as indicated in the figure legend, scaled by factors of 10 
between each defect shape to avoid nearly overlapping in the 
figure. The long-dashed lines are a guide for the eye along the 
numerical results. The TFAs results of Eq. (|14|) are the solid 
lines. The square-well results are compared with the analyt- 
ical expression of Carr et ai, Eq. (|15|) . as the short dashed 
line. 



of the square-well trap, we also show the approximate 
analytical solution given by Eq. (fT5|) (short-dashed line) . 

The reason to plot gi,maxRo "vs. i?oV^ in a log-log 
scale and not gi,maxRo/ /{RoV^) is to show the range 
of values that gi^max might take over the range we used 
in Ro^/\^, and not just to have a constant value defined 
by the shape factor. 

As observed in Fig. [21 some symbols are seen to over- 
lap, confirming the basic scaling rule for different widths 
and depths of the same defect. For large Ro\/Vo the 
numerical results follow closely the TFA. This is due to 
large Ro\/Vo which requires a large number of atoms to 
fill up the defect, i.e., the system is in a strong nonlin- 
earity regime where the TFA is valid. Note though, that 
this also means that our approximation of holding the 
wavefunction in a tight transverse state will eventually 
breakdown. For values of Ro\/Vq < 1, the kinetic energy 
term starts to dominate since gi^max gets smaller. It is in 
this region that the difference between the defect shapes 
is revealed. However, this is also the region where the 1- 
D NLSE itself is not valid anymore as mentioned in the 
introduction. 

As the only difference between the 1-D defects is their 
shape factors [Eq. (fTT])]. we can also plot all of them in 
a single figure. Figure |3] shows the results by plotting 
gi,maxRo/oiiD as a function of Rt^^/Va. For _Ro\/Vo > 1 
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FIG. 3. Scaling laws for gi.maxRo/aiD as a function of RoVVo 
for the five different 1-D defect shapes. The solid line is the 
TFA which becomes identical for all defects. The symbols are 
the numeric results as obtained from Eq. (|16[l. The short- 
dashed line is the approximation of Carr et al. [2^ for the 
square- well. 



there is barely any difference for all of the trapping de- 
fects in the factor gi^maxRa/aiD, confirming our gener- 
alized scaling rule. Thus, under these conditions, the 
square- well defect holds the most atoms. It is followed 
by the circular defect, then by the harmonic well defect 
and the double harmonic well, and finally the triangular 
defect will hold the least atoms for the same width and 
depth defects. This is in both the numerical results and 
the TFA [Eq. p^]. 

In Fig. 01 we show the scaled wavefunctions, ^(x) = 
^/Ro'i'{x), for the square- well defects for two values of 
gi,maxRo, and two values of RoV^^ as a function of a; = 
x/Rq. The two cases shown give us a small and a large 
value of gi^maxRo- For the first case, RoV^o = 1.5811 
and gi^maxRo = 6.284, in particular, the two square-well 
defect wavefunctions have Rf) = 5.0 and Vq = 0.1 o.u. 
(solid line) for one case, and _Ro ~ 0.5 and Vq — 10.0 
o.u. (* symbol) for the other case. That is, a wide and 
shallow potential trap vs. a tight and deep potential, but 
both with small Ra^/\^ value. 

We note in Fig. 21 that the scaled wavefunction show 
the same tunneling for both cases of Rq\/Vo. However, 
due to the factor ^/Rq in front of ^I*, both cases would 
have different tunneling in the x oscillator units space. 
The second case considered in Fig. 21 shows the wavefunc- 
tion of two square-well defects with Rq^/\^ ~ 14.1421 
and gi.maxRa — 415.770, in particular, Rq = 20.0 and 
Vo = 0.5 o.u. (dashed line) and Rq = 2.0 and Vq = 50.0 
o.u. (• symbol). Again, we have a wide and shallow po- 
tential trap vs. a tight and deep trap, but now with a 
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FIG. 4. (Color on-line). Scaled 1-D numerical wavefunctions, 
^(x) = V Ro'^{x/Ro), for various square-well defects. The 
solid line corresponds to Ro = 5.0 and Vo ~ 0-1 o.u.; the 
blue (*) symbols correspond to Ro = 0.5 and Vo = 10.0 o.u., 
such that both of them have the same Ro\/Vo = 1.5811 and 
gi,maxRo = 6.284. The dashed line denotes _Ro = 20.0 and 
Vo = 0.5 o.u., while the (•) symbol is for Ro = 2.0 and Vo = 
50.0 O.U., both of these wavefunctions have Ro^/Vo ~ 14.1421 
and gi^maxRo = 415.770, thus, showing the same reduced 
wavef unction. 



large values for i?o\/ Vq. In this case we are in the high 
i?o\/Vo^ region where we are closer to the TFA as shown 
by the shape of the wavefunction which is closer to the 
potential shape. 

The scaUng was also verified for the other potentials. 
In Fig. [51 we show the corresponding scaled wavefunc- 
tions for the circular (a), harmonic (b), double harmonic 
(c) and triangle (d) trap potentials for the same width 
and depth parameters as those shown in Fig. |4l Note 
once again how the solutions to the NLSE satisfy the 
scaling rule, showing the same reduced wavefunction in 
the reduced units. For the case of the double harmonic 
trap and the triangle potential trap [Fig. ((St) and ([5jl)] 
the wavefunction smooths out for small gi^maxRo values 
in the regions where the potential is non-differentiable in 
contrast to the strong interaction region where the wave- 
function shows the same shape as the trapping potential. 
Thus, we have confirmed for a given defect that two dif- 
ferent shapes with the same Ro\/Vo will have the same 
gi,maxRo and the same reduced wavefunction. 
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FIG. 5. (Color on-line). Scaled 1-D numerical wavefunctions, 
$(a;), for the (a) circular, (b) harmonic, (c) double harmonic 
and (d) triangle potential. The labels and parameters are the 
same as in Fig. [D 



we show the results for Rx — Ry = Rq ~ 1.0 o.u. and 
Vo — 50.0 o.u. for the results in red lines (solid lines). 
Interestingly, the wavefunction takes a shape similar to 
the potential the particles are being held, in the same way 
as the 1-D case. In the same figure, we show the results 
for Rj. = Ry — Ro — 5.0 o.u. and Vq — 2.0 o.u. for 
the results shown by the symbols (blue squares). Both 
cases have Roy/Vo = 7.071, therefore both of them have 
the same scaled wavefunction, $ = y/R^Ry'^f = Rq^, 
confirming our scaling rule, i.e. the same reduced density 
profile. 

In Fig. [7j we show the ratio of the nonlinear cou- 
pling term g2,max to the shape form factor a2D for the 
2-D results as a function of Ro^/Vo. The solid lines 



are the Thomas-Fermi results and the symbols are the 
data obtained by our numerical procedure. Again, for 
RoVVq > 10 the Thomas- Fermi results follow closely the 
numerical data for g2,max/oi2D showing a universal be- 
havior for the non-linear coupling term. Discrepancies 
start to appear in the results for small values of i?o\/Vo, 
dependent of the defect shape. This is a consequence of 
the neglect of the kinetic energy term in the TFA, and 
thus the neglect of tunneling for these weakly bound par- 
ticles. 



D. The 3-D trap cases 



C The 2 D trap cases '^'^ order to visualize the density profile of a 3-D wave- 

function, we project the reduced wavefunction onto the 

„ , ,. . , 1 • r?n 1 z-axis such that 

For the two-dimensional cases we show m Fig. |6j the 

probability density for the three potential cases (square, _ 2 /" i /- - -m2j- 

harmonic, and pyramid well potentials). In these cases \^\.x^y)\ ^ J \^[x,y, z)\ dz . (43) 
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FIG. 6. (Color on-line). Numerical scaled density profile solu- 
tions, |<l>(x, y)p for the 2-D NLSE as a function of x and y for 
the square, parabolic, and pyramid trap defects for Vo — 50.0 
and Rx = Ry = Ro = 1.0 o.u. for the solid lines (red lines). 
In the same figure we show the results for Vo ~ 2.0 and 
Rx = Ry = Ro ~ 5.0 o.u. represented by the squared symbols 
(blue), thus confirming the scaling density profile rule for the 
2-D case. 



Thus, for the three-dimensional case, we show in Fig. 
[5] the projected scaled density profiles, |<f>(a;,y)p, for 
the spherical, rectangular cuboid and harmonic well po- 
tentials. For these two cases we show the results for 
= Ry = Ro = 1.0 o.u. and Vq = 50.0 o.u.. The 
wavefunction takes a shape similar to the defect that the 
particles are being held, as the 1- and 2-D wavefunctions 
also tended to do. In the same figure, we show the results 
for Rx — Ry — Rq = 5.0 o.u. and Vq — 2.0 o.u.. Both 
cases have Roy^ = 7.071, and therefore have the same 

I 3/2 

scaled wavefunction $ — R^RyR^^ = R^ ^f, as well 
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FIG. 7. Scaling law for the 2-D defects, g2,max / 0120 of the 
nonlinear coupling constant as a function of Rq\/Vo for the 
square, parabolic, and pyramid trapping potential for the 2- 
D NLSE. The solid line is the Thomas-Fermi result. Symbols 
are our numerical results at the threshold of delocalization. 
The upper/middle/lower symbols are, respectively, for square, 
parabolic and pyramid shaped defects. 



as the same gs.max/ Ro confirming that our scaling rule 
further holds in 3-D. 

In Fig. |9l we show the ratio of the nonlinear coupling 
term ga^max to the shape form factor a^o for the three 
dimensional results as a function of RQ^/^). The solid 
straight line is the Thomas-Fermi results and together 
with the symbols are the data obtained by our numerical 
procedure. Again, for Rq\/\^ > 10 the Thomas- Fermi re- 
sults follow closely the numerical data for g^,max/ ct^oRo 
showing an universal behavior for the non-linear coupling 
constant. For small values of Rq^/Vq, discrepancies, de- 
pendent of the potential trap shape, start to appear in 
the results consequence of the kinetic energy term. In 
the same figure, we show the variational results for Ad- 
hikari [24] . As mentioned in the introduction, his results 
overestimates the TFA or our numerical results. 



E. Application: A 3-D cube trap 

The application and implications of this work follows 
directly from the scaling law for the non-linear coupling 
constant and the reduced density profile for the 1-, 2-, 
and 3-D cases considered. The NLSE in physical units is 

i 

{-|^V2 + y,.,,(r) + (44) 

^ ^-0 r ' \ ^-0 r = £0*0 r . 

m J 




FIG. 8. (Color on-line). Numerical projected scaled density 




FIG. 9. Scaling law of the 3-D nonlinear coupling constant 
divided by the potential shape factor a^o as a function of 
RoV^ for the spherical, rectangular cuboid, and parabolic 
trapping potential for the 3-D NLSE. The solid line is the 
Thomas- Fermi result. The symbols are our numerical results 
near delocalization. The upper/middle/lower symbols are, 
respectively, for spherical, rectangular cuboid, and parabolic 
shaped defects. The dashed line represents the variational 
results from Adhikari '53] for the spherical well. 



where 



53 



47ra,(A^- 1) 



V = 



m(l 



xlyl 



)l/3 
)2/3 



(46) 



(47) 
(48) 



As an example we again choose ^^Rb atoms with = 
100 flo as in the experiment of Ref. [H. We consider the 



profile solutions I $(x y) I to the 3-D NLSE as a function of ^^^^ ^ ^^^^ ^.^^^ ^^^^^^ ^^^j^ I. = ly = h = 



and y for the spherical, rectangular cuboid, and parabolic well 
trap potentials for = Ry = Rz = 1.0 o.u. and Vo — 50.0 
for the solid lines (red lines). In the same figure we show the 
results for Rx = Ry = Rz ~ 2.0 o.u. and Vo ~ 5.0 by the 
squared symbols (blue), thus confirming the scaling density 
profile rule for the 3-D case. 



For a trap with characteristic harmonic oscillator lengths 
Ix, ly, and Iz for each Cartesian coordinate, such 
that a transformation to oscillator units requires r = 
{IJyhflH, and the NLSE in o.u. becomes 



1.222 X IQ-'^ cm = 23000ao and thus a volume of 8 x Z^. 
The TFA then gives us = 0.05464(7V - 1) and V = 
2.68 X 10* Vtrap{o-n/ K), where the potential trap depth 
is given in units of absolute temperature (Kelvin). Thus, 
for i?o covering the range from 1 to 10 o.u., as used in Sec. 
lIVDI with physical values that correspond to a BEC with 
extension between 1.22 and 12.2 /im. Similarly, for V — 
Vb covering the range from 1 to 50 o.u. this is equivalent 
to traps depths of 3 to 150 nK, which are values well 
within the range of the experiment. The values reported 
in Fig [3] thus correspond to a range from 200 to 2 x 10^ 
atoms, which are typically found in experiments. 



CONCLUSIONS 



--^' + Vtrap{v)+9:i\^{v)Y 



^{v) = e^ir) (45) 



In conclusion, by means of the Thomas-Fermi approx- 
imation at the delocalization threshold, ei = 0, one can 



easily determine an approximation to the maximum non- 
linear interaction strength Qi^max of the Gross-Pitaevskii 
equation and therefore, the maximum number of atoms, 
Nmax, that can be trapped by i =1-, 2-, and 3-D poten- 
tials that contain a defect. 

The scahng laws of Eqs. dH]), and ([35]) are de- 

pendent on both the depth and length of the potential, 
and remain valid over a surprisingly vifide-range of pa- 
rameters. The TFA relies on a vifavefunction that is con- 
strained by the shape of the potential, never diffusing 
into the classically forbidden region. At the Qi^max point 
the TF wavefunction becomes unbound, and that this 
mimics the actual solution to the NLSE which diffuses 
towards infinity as gi approaches gi^max is remarkable. 
It would, however, be worthwhile to extend this work 
through a more sophisticated variational treatment [33| . 

The existence of these scaling laws is useful because it 
captures in a simple expression the number of atoms of a 
BEC that can be trapped by a potential defect on waveg- 
uide or in free-space. It would be interesting to further 
examine the dynamic scattering /trapping of a BEC as 
it propagates through such defects [1^, [20] • Essentially, 
the non-linear term enables a continuum of non-linear 
bound states, as opposed to the quantized single-particle 



eigenstates. For example, if N^ax = 200 atoms, then any 
number lower than that will also be able to be trapped. 
How easy is it for a BEC with N > Nmax to fill up the 
defect as it goes past remains to be demonstrated, and 
our initial calculations within the NLSE show very little, 
if any, filling up of the defect. It will also be interest- 
ing to consider that the presence of any localized atoms 
will tend to smooth out the defect potential that a con- 
sequent BEC interacting with the defect will experience. 
This may have an impact, e.g. on Anderson-type local- 
ization. 
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